function [] = readProviderMovements(inFolder, outFolder)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% This script takes as input the following file:
% inFolder/providerData.csv:        UserId  | providerCreationDate_julian | providerDeletionDate_julian | CantonId
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Add paths
addpath('../code/general_funcs');

%%%%% Define geographic levels
GEO_VARS = {'CantonId'};
NumGeoVars = length(GEO_VARS);

%%%%% Define paths and input files
file1 = sprintf('%s/providerData.csv', inFolder);
geoDefFile     = sprintf('%s/../covariates/Scope_CantonId.csv', inFolder);
dayLabelsFile  = sprintf('%s/../covariates/periodLabels_day.csv', inFolder);
weekLabelsFile = sprintf('%s/../covariates/periodLabels_week.csv', inFolder);

%%%%% Read input files
% Read geo-def file and get geo codes that will allow me to convert from geoCode to integer
opts_geo = detectImportOptions(geoDefFile);
for GEO_VAR = GEO_VARS;
	opts_geo = setvartype(opts_geo, GEO_VAR, 'char');
end
Mgeo = readtable(geoDefFile, opts_geo);

%% Get period labels (days and weeks)
M1 = readtable(dayLabelsFile);
dayLabels = table2cell(M1(:,2));
dayLabels2 = cellfun(@(x) datestr(x, 'yyyy-mm-dd') , dayLabels,'UniformOutput',false);
clear M1;

M2 = readtable(weekLabelsFile);
weekLabels = table2cell(M2(:,2));
clear M2;

% Read providerData file
opts = detectImportOptions(file1);
for GEO_VAR = GEO_VARS;
	opts = setvartype(opts, GEO_VAR, 'char');
end
M3 = readtable(file1, opts);

%%%%% Get all labels: providerLabels, geoCodes
% Make providerLabels and NumProviders
providerLabels  = sort(make_integers(table2array(M3(:,1)))); % NumProviders x 1
NumProviders = length(providerLabels);

% Make geoCodes
geoCodes = cell(1,NumGeoVars);
for gidx = 1:NumGeoVars
	GEO_VAR = GEO_VARS{gidx};
	[~,mycolidx] = ismember(GEO_VAR, Mgeo.Properties.VariableNames);
	geoCodes{gidx} = unique(table2cell(Mgeo(:,mycolidx)));
	clear mycolidx;
end

% Get car2geo, user2geo, provider2geo
provider2geo  = a_9_getGeoMapping(M3, GEO_VARS, geoCodes, 'UserId', providerLabels);  % cell(1,NumGeoVars) of NumProviders x 1  --> gives geo index of each provider

%%%%% Get provider2creationDate and provider2deletionDate: NumProviders x 1 --> gives real number representing timestamp
provider2creationDate  = a_9_getNumericVariableFromTable(M3, 'UserId', providerLabels, 'providerCreationDate_julian');
provider2deletionDate  = a_9_getNumericVariableFromTable(M3, 'UserId', providerLabels, 'providerDeletionDate_julian');

%%% Clean up
clear file1 opts opts_geo;
clear M3 Mgeo ;

%%%%% Put together carCreations, carDeletions, userCreations by location and time

% Loop over time aggregation levels
timeAggregationLevelNames = {'week', 'day'};
for tidx = 1:length(timeAggregationLevelNames)
	
	timeAggLevel = timeAggregationLevelNames{tidx};
	if strcmp(timeAggLevel, 'day')
		Tmax = 1826;
		timeLabels = dayLabels;
		timestamp_2_tt_function = @(time_val) ceil(time_val);  % Discretizes timestamp into tt_vec (time stamp to discrete day idx)
	end;
	if strcmp(timeAggLevel, 'week')
		Tmax = 260;
		timeLabels = weekLabels;
		timestamp_2_tt_function = @(time_val) ceil(time_val/7); % Discretizes timestamp into tt_vec (time stamp to discrete week idx)
	end;

	for gidx = 1:NumGeoVars
		GEO_VAR = GEO_VARS{gidx};
		geoLabels = geoCodes{gidx};
		K = length(geoLabels);
		provider_kk_vec  = provider2geo{gidx};  % NumProviders x 1  --> gives index between 1 and K
		
		outputFile = sprintf('%s/providerMovementsData_%s_%s.mat', outFolder, GEO_VAR, timeAggLevel);
		
		disp(sprintf('Saving file: %s...', outputFile));
		
		%%%%% Count events by time and location (and past events by location)
		% Provider creations
		tt_vec = timestamp_2_tt_function(provider2creationDate); % NumProviders x 1
		[pastProviderCreations, providerCreations]   = a_9_countEvents_ByTimeAndLocation(tt_vec, provider_kk_vec, Tmax, K);
		
		% Provider deletions
		tt_vec = timestamp_2_tt_function(provider2deletionDate); % NumProviders x 1
		[pastProviderDeletions, providerDeletions]   = a_9_countEvents_ByTimeAndLocation(tt_vec, provider_kk_vec, Tmax, K);
		
		%%%%% Get initial base
		initialBase_providers  = pastProviderCreations - pastProviderDeletions; % 1 x K
		
		%%%%% Get net creations
		providerNetCreations = providerCreations - providerDeletions; % Tmax x K
		
		%%%%% Get base at END of period tt
		base_providers  = initialBase_providers + cumsum(providerNetCreations); % Tmax x K
		
		%%%%% Get base at START of period tt
		laggedBase_providers  = base_providers  - providerNetCreations; % Tmax x K
		
		%%%%% Combine everything together
		providerCreations     = reshape(providerCreations,   [Tmax*K 1]); % (Tmax*K) x 1
		providerDeletions     = reshape(providerDeletions,   [Tmax*K 1]); % (Tmax*K) x 1
		laggedBase_providers  = reshape(laggedBase_providers,  [Tmax*K 1]); % (Tmax*K) x 1
		
		movementVarNames = {'NewProviders', 'DeletedProviders', 'LaggedProviders'};
		movementVars = [providerCreations providerDeletions laggedBase_providers]; % (Tmax*K) x NumVars
		movementVars = reshape(movementVars, [Tmax K size(movementVars,2)]); % Tmax x K x NumVars
		
		movementData.values        = movementVars; % Tmax x K x NumVarZ
		movementData.VarNames      = movementVarNames;
		movementData.geoLabels     = geoLabels;
		movementData.timeLabels  = timeLabels;
		
		%%%% Output to file
		save(outputFile, 'movementData', '-v7.3');
	end
end


%%%%% Output continuous-time data
for geoidx = 1:NumGeoVars
	GEO_VAR = GEO_VARS{geoidx};
	geoLabels = geoCodes{geoidx};
	providerCreations_ctsData  = [provider2creationDate  provider2geo{geoidx}  1*ones(NumProviders,1)]; % NumProviders x 3
	providerDeletions_ctsData  = [provider2deletionDate  provider2geo{geoidx}  2*ones(NumProviders,1)]; % NumProviders x 3
	providerDeletions_ctsData  = providerDeletions_ctsData(~isnan(provider2deletionDate),:);  % NumProvidersDeleted x 3
	
	eventsData = [providerCreations_ctsData; providerDeletions_ctsData]; % NumEvents x 3
		% Column 1: time stamp (julian day with ref. of 2013-01-01)
		% Column 2: kk (location)
		% Column 3: event type: (1=provider creation, 2=provider deletion)
	
	% Sort eventsData by increasing timestamp
	[~,sortidxes] = sort(eventsData(:,1));
	eventsData = eventsData(sortidxes,:);
	
	timestampsData.geoLabels = geoLabels;        % K x 1
	timestampsData.timestamps = eventsData(:,1); % NumEvents x 1
	timestampsData.kk_idxes   = make_integers(eventsData(:,2)); % NumEvents x 1 (values between 1 and K)
	timestampsData.jj_idxes   = make_integers(eventsData(:,3)); % NumEvents x 1 (code between 1 and NumEventTypes, that indicates even type)
	timestampsData.jj_labels = {'NewProviders', 'DeletedProviders'}; % 1 x NumEventTypes
	
	outputFile = sprintf('%s/providerMovementsData_%s_ctsTime.mat', outFolder, GEO_VAR);
	disp(sprintf('Saving file: %s...', outputFile));
	save(outputFile, 'timestampsData', '-v7.3');
end

end